Learning Objectives
- Introduce map outlines available in the
mapspackage.
- Learn to set and transform coordinate systems using the
sfpackage.
- Reshape and plot spatial data using the
tidyverse.
- Plot raster data using the
terrainrandrasterpackages.
- Add scale bars and north arrows using the
ggspatialpackage.
- Next week: Apply skills to create a map of field sites or sampling locations.
Today we will cover the basics of creating maps in R. By the end of
this workshop, you will be able to create a map and add labeled points,
raster data, and shape files. Our introduction to mapping will, in many
ways, be an extension of data manipulation and plotting within the
tidyverse.
Open your Environmental Data Analysis project, and create a new script named “A10_Mapping”. Download the mapfiles folder from Brightspace, and extract the files into a new folder named “mapfiles” in your project folder. Remember to back up your project folder after completing the workshop today.
Load the packages you will need today. Computers 1-20 in Hudson 032 should have all of these packages installed already. If you are working on your personal computers, you will need to install the new mapping packages.
#Load packages
library(tidyverse)
library(maps)
library(sf)
library(terrainr)
library(raster)
library(ggspatial)
The first thing we will do is download a map of the continental
United States. The maps package contains a file named
“state” that contains georeferenced polygons for all of the lower 48
states. We will use the map() function to import these data
and assign them to the object state. The
fill = T argument specifies that you would like each state
to be stored as a separate polygon rather than a series of lines.
# Get states
state <- map("state", fill = T)
This function should return an outline of the continental United
States. Note: The package purrr (part of the
tidyverse) and the map package both contain
the function map(). This may cause an error when you
attempt to retrieve the map data. When you are working with packages
that contain functions of the same name, you can specify which package
you mean by adding the package name ahead of the function using the
syntax package::function(). So for the example above, you
would specify that you want to use the map function from the maps
package using maps::map("state", fill = T).
If you look at your state object in the environment
window, you will see that the function has created a list with four
vectors. These vectors contain x and y coordinates (longitude and
latitude), the bounds of the area, and the names of the states. We would
like to turn this list into a spatial object that we
can use to visualize and analyze relationships between locations. The
sf package allows us to do just that. We will use the
st_as_sf() to transform our state map object
into a spatial data frame.
# Create sf object
statesf <- st_as_sf(state)
Now if you take a look at the statesf object in the
environment window, you will see that the function has created a data
frame with each state as a separate row in the ID column
(including the District of Columbia) and a geom column that
contains a list specifying the x-y coordinates of the outline of the
state polygon. Run your object name to see more information that is
stored in this object.
statesf
## Simple feature collection with 49 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -124.6813 ymin: 25.12993 xmax: -67.00742 ymax: 49.38323
## Geodetic CRS: +proj=longlat +ellps=clrk66 +no_defs +type=crs
## First 10 features:
## ID geom
## alabama alabama MULTIPOLYGON (((-87.46201 3...
## arizona arizona MULTIPOLYGON (((-114.6374 3...
## arkansas arkansas MULTIPOLYGON (((-94.05103 3...
## california california MULTIPOLYGON (((-120.006 42...
## colorado colorado MULTIPOLYGON (((-102.0552 4...
## connecticut connecticut MULTIPOLYGON (((-73.49902 4...
## delaware delaware MULTIPOLYGON (((-75.80231 3...
## district of columbia district of columbia MULTIPOLYGON (((-77.13731 3...
## florida florida MULTIPOLYGON (((-85.01548 3...
## georgia georgia MULTIPOLYGON (((-80.89018 3...
Note that the range from the maps object appears as the
“Bounding box” in the sf object, the geometry appears as a
MULTIPOLYGON, and the coordinate system of the map appears
as the Geodetic CRS. In this case, the coordinate reference
system (CRS) for our map is WGS 84, which is commonly
used in North America.
Once you have created an sf object, you can use
tidyverse functions to reshape and plot the data. Let’s use
the ggplot() function to create a quick map. We will use
ggplot() to create the plot object and identify our data
frame, and then add a geom_sf() layer to add our sf object
geometry to the map.
ggplot(statesf) +
geom_sf()
Note You may notice that your degree symbols are not plotting correctly. If so, this is a graphics issue. You can fix this issue in RStudio by selecting
Tools->Global Options. This will open a window. Choose theGraphicstab at the top of the window, and change theBackendto “Cairo.”
All of the arguments in ggplot can be added to customize your plot.
You can change the color of the line and the fill of the polygons, or
the width of the outline. You can change the theme using the
theme_minimal() or the theme_bw() functions we
have already learned. A common theme for mapping is
theme_void, which removes all elements except for the map
layers.
ggplot(statesf) +
geom_sf(color = "#283044", fill = "#74D3AE", size = 0.5) +
theme_void()
Note that fixed modifiers are included outside of a
aes() mapping function. If you wished to add a modifier
based on a variable, that would need to be included inside the
aes() function. For example, if I wanted to plot every
state as a different color, I could add the argument
aes(fill = ID) to the geom_sf() layer.
ggplot(statesf) +
geom_sf(aes(fill = ID)) +
theme_void()
You can also easily plot points on the map based on coordinates in a
tibble or data frame. When you add a new layer, in this case a
geom_point() and a geom_text() layer, you need
only specify the data frame where the information is located using
data = and then specify your aesthetics. Let’s add points
to show the positions of SUNY Plattsburgh and Washington D.C. on our
map.
First we will create a tibble with names and coordinates for the two locations, then we will add points and labels to our map.
# Create table with coordinates and names
labs <- tibble(
long = c(-73.464279, -77.015541),
lat = c(44.692614, 38.892753),
names = c("SUNY Plattsburgh", "Washington D.C."))
# Add point and text layers to map
ggplot(statesf) +
geom_sf(color = "#F5DBCB", fill = "#74D3AE", size = 0.5) +
geom_point(data = labs, aes(x = long, y = lat),
color = "#283044", size = 5) +
geom_text(data = labs, aes(x = long, y = lat, label = names),
hjust = 1, nudge_x = -1, color = "#283044")
Notice that I added the arguments hjust and
nudge_x to orient the labels so that they are readable. If
you are able to install packages, the package ggrepel has a
function geom_label_repel() that you can add in place of
geom_text(). This function will automatically space labels
so that they are readable.
Challenge 1
Create a plot of the United States with every state plotted as a different color.
Hint: Look at thestatesfobject to determine which variable to use. Think about where you need to add thefill =argument to change fill based on a variable in thestatesfdata frame.
Change the color that is used to outline the states.
Consult the ggplot2 cheat sheet or Google to figure out how to exclude the legend from appearing on your map.
When you are satisfied with your map, useggsave("A10_C1.jpg")to export your map image.
Upload your plot to Challenge 1 on Brightspace.
ggplot(statesf) +
geom_sf(aes(fill = ID), color = "white") +
theme_minimal() +
theme(legend.position = "none")
ggsave("figures/A10_C1.jpg")
## Saving 7 x 5 in image
Because the sf object is stored as a data frame, we can
use the tidy tools of dplyr to retain just certain parts of
it. For example, we can select just NY and VT and then plot them.
# Create data frame for just NY and VT
ny_vt <- statesf %>%
filter(ID %in% c("new york", "vermont"))
#Plot only those two states
ggplot(ny_vt) +
geom_sf(fill = "#74D3AE", color = "#F5DBCB") +
theme_minimal()
So far, we’ve added point and text data to our base map by adding
layers to our ggplot and specifying a data object. We can use a similar
approach to layer additional spatial data onto our map. In the next
section, we will import a shape file for Lake Champlain as an
sf object, transform it into the same coordinate reference
system (CRS) as our current layers, and then plot it using the
geom_sf() function.
A shape file is a data storage format that stores the location, shape, and various attributes of spatial features. A shape file is stored as a set of related files (i.e., cpg, dbf, prj, shp, and shx) with the same file name prefix. All of these files are needed to plot the geometry appropriately and should be stored in the same folder in your file directory.
Let’s import our shape file for Lake Champlain from the
mapfiles folder using the function
st_read() in the sf package. The syntax for
this function is st_read("folder name", "file prefix").
# Import shape file data from mapfiles folder
lc <- st_read("mapfiles", "VT_Lake_Champlain")
## Reading layer `VT_Lake_Champlain' from data source
## `C:\Users\malld001\Dropbox\PlattsburghCourses\ENV422_EnvDataAnalysis_Spring2024\A10Mapping\mapfiles'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 9 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 420142.5 ymin: 115042.3 xmax: 454613.1 ymax: 287432.3
## Projected CRS: NAD83 / Vermont
Take a look at the attributes for the sf object you have
just imported. Note that the Projected CRS differs from the
“WGS 84” system on the maps we have made up to this point. We should
transform the coordinate system of our shape file if we want the layers
to match up properly.
When we plot a map in two dimensions, we are necessarily translating 3D information into a 2D picture. This is easier for small areas than large areas, but maps are always a projection of a 3D reality. R can read more than 6,000 different projections and transform them to a different CRS. Yes there are that many! Fortunately, within a region, only ~3-5 are typically used. Most spatial data are stored with this information, or the information appears in the meta-data. For example, our Lake Champlain shape file specifies the CRS as NAD83 / Vermont.
In some cases, you may import (or be given) a file with no CRS associated with it. You will need to figure out the projection system used in order to transform it to the projection you are using. Some of the most commonly used systems in North America are listed in the table below.| CRS | EPSG | Note |
|---|---|---|
| WGS84 | 4326 | Commonly used by organizations that provide GIS data for the entire globe or many countries. CFS used by Google Earth. |
| NAD83 | 4269 | Most commonly used by US federal agencies |
| NAD27 | 4267 | Old version of NAD83 |
| Mercator | 3857 | Tiles from Google Maps, Open Street Maps, Stamen Maps |
| UTM, Zone 10 | 32610 | See map below |
| UTM, Zone 11 | 32611 | See map below |
| UTM, Zone 12 | 32612 | See map below |
| UTM, Zone 13 | 32613 | See map below |
| UTM, Zone 14 | 32614 | See map below |
| UTM, Zone 15 | 32615 | See map below |
| UTM, Zone 16 | 32616 | See map below |
| UTM, Zone 17 | 32617 | See map below |
| UTM, Zone 18 | 32618 | See map below |
| UTM, Zone 19 | 32619 | See map below |
By Chrismurf at English Wikipedia, CC BY 3.0, https://commons.wikimedia.org/w/index.php?curid=40690482
To transform our Lake Champlain shape file from NAD83 to WGS84, we
will use the st_transform() function in the sf
package. We will name the transformed shape file object
lcsf.
lcsf<- lc %>%
st_transform("WGS84")
Run the lcsf object and plot it to confirm that the
transformation was successful and that the geometry of the shape looks
correct.
# Print lcsf shape file object
lcsf
## Simple feature collection with 1 feature and 9 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -73.48796 ymin: 43.53135 xmax: -73.07608 ymax: 45.08523
## Geodetic CRS: WGS 84
## OBJECTID COMID GNIS_NAME AREASQKM REACHCODE FTYPE GDB_GEOMAT
## 1 29848 <NA> Narrows, The 1146.091 02010008009170 LakePond <NA>
## SHAPESTAre SHAPESTLen geometry
## 1 1141293318 1102278 POLYGON ((-73.14224 45.0851...
# Plot lcsf
ggplot(lcsf) +
geom_sf()
Now that our shape file is in the right coordinate system, we can add
it to our map of New York and Vermont by adding an additional
geom_sf() layer.
ggplot(ny_vt) +
# add state outlines
geom_sf(fill = "#74D3AE", color = "#F5DBCB") +
# add lake champlain shape
geom_sf(data = lcsf, color = "#1282A2", fill = "#1282A2") +
# remove background
theme_void()
So transforming the Lake Champlain shape file to the correct coordinate system was fairly straightforward. Let’s take a look at a trickier case. We’ll import the shape file that represents the extent of the 2018 forest fire at the Altona Flat Rock pine barrens.
First, we’ll import the shapefile from the mapfiles folder using the
st_read() function.
# Import shape file data from mapfiles folder
fr <- st_read("mapfiles", "fire_perimeter")
## Reading layer `fire_perimeter' from data source
## `C:\Users\malld001\Dropbox\PlattsburghCourses\ENV422_EnvDataAnalysis_Spring2024\A10Mapping\mapfiles'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -8198195 ymin: 5599219 xmax: -8195719 ymax: 5602464
## CRS: NA
Note that no coordinate system is listed under CRS. Moreover, the
coordinates for the bounding box look very different from the latitudes
and longitudes we’ve seen so far. That means we’re dealing with a very
different coordinate system, and we don’t know what it is! We will need
to use the st_set_crs() function to specify the current
coordinate system before we can transform the data to the “WGS84” system
that we are using. The syntax for this function is
dataframe %>% st_set_crs(EPSG)
After much trial and error (set -> transform -> plot ->
work?), I determined that the coordinates were in Mercator (EPSG =
3857). We will use st_set_crs() to set the object to its
current projection and then st_transform() transform it to
WGS84. We’ll save the transformed shape file as frsf.
frsf <- fr %>%
st_set_crs(3857) %>%
st_transform("WGS84")
frsf
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -73.64564 ymin: 44.85816 xmax: -73.62339 ymax: 44.87882
## Geodetic CRS: WGS 84
## Id geometry
## 1 0 POLYGON ((-73.64564 44.8639...
Note that the new frsf object has a coordinate system
assigned, and that the Bounding box now looks like it contains latitude
and longitudes in a sensible region of the planet. Try plotting the
shape file on its own and on the ny_vt map to see how it looks at the
different scales!
You can also combine layers available in the maps
package by adding geom_sf() layers just as we did for the
shape file.
In this section, we will create a map of New York with the boundaries
of counties included. First we will use map() to import the
county polygons and st_as_sf() to transform the data to a
spatial data frame, just as we did for states.
# Get counties from maps package
county <- maps::map("county", fill = T, plot = F)
# Create sf object
countysf <- st_as_sf(county)
Take a look at the countysf object. Notice that the ID
contains a single column with both the county and state names. As we
learned in the data management section, this format is not ideal.
However, we can use functions in the dplyr package to
easily separate this column into two columns, one for state and one for
county.
# Overwrite data frame with new one that contains separate state and county columns
countysf <- countysf %>%
# create a new column, ID2, that is identical to ID
mutate(ID2 = ID) %>%
# separate the new column into two columns, using "," as the delimiter
separate(ID2, into = c("state", "county"), sep = ",")
When you view the countysf data frame now, you will see
the original column as well as your new separated columns containing the
state and county names.
Now let’s use the tidyverse to select only counties in
New York and create a map.
countysf %>%
# select only counties in new york
filter(state %in% c("new york")) %>%
# create plot and add geom_sf layer
ggplot() +
geom_sf(fill = "#74D3AE", color = "#F5DBCB") +
# eliminate background features
theme_void()
Once we have this map, we can add as many layers as we like. As an example, let’s add a layer with the outline of New York.
# Create a data frame with just the outline of New York from the statesf data frame
ny <- statesf %>%
filter(ID %in% "new york")
# Create a data frame that contains only the county boundaries in NY from the countysf data frame
county_ny <- countysf %>%
filter(state %in% c("new york"))
# Copy map code from above
ggplot(county_ny) +
geom_sf(fill = "#74D3AE", color = "#F5DBCB") +
theme_void() +
# add layer for outline of New York state
geom_sf(data = ny, color = "#283044", fill = NA)
Challenge 2
Based on what you have learned, create a map for the state of New York that includes county boundaries and a polygon showing the outline of Lake Champlain.
Add a layer that places a labelled point over SUNY Plattsburgh.
When you are satisfied with your map, useggsave("A10_C2.jpg")to export your map image.
Upload your plot to Challenge 2 on Brightspace.
# Create an object that contains only the county boundaries in NY
county_ny <- countysf %>%
filter(state %in% c("new york"))
# Create map
ggplot(county_ny) +
# add county boundaries
geom_sf(fill = "#74D3AE", color = "#F5DBCB") +
theme_void() +
# add layer for outline of New York state
geom_sf(data = ny, color = "#283044", fill = NA) +
# add layer for outline of Lake Champlain
geom_sf(data = lcsf, color = "#1282A2", fill = "#1282A2") +
# add point at SUNY Plattsburgh coordinates
geom_point(aes(x = -73.464279, y = 44.692614), size = 3, color = "#283044") +
# add text at SUNY Plattsburgh coordinates and adjust position
geom_text(aes(x = -73.464279, y = 44.692614, label = "SUNY Plattsburgh"),
hjust = 1, nudge_x = -0.2, color = "#283044") +
# remove background elements
theme_void()
ggsave("figures/A10_C2.jpg")
## Saving 7 x 5 in image
Now that we can create a nice map of the counties, ideally we would
also like to add data. To do this, we will first import a data frame
that I downloaded from Wikipedia and cleaned in Excel. We will import
this file from the mapfiles folder that you downloaded from
Brightspace. When you view the data frame, note that I added a column
for “county” that has the name of the counties exactly as they appear in
the county_ny data frame.
#Import wikipedia data
ny_wiki <- read_csv("mapfiles/nycounty_info.csv")
## Rows: 62 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): county, CountyName, CountySeat, formed_from, named_for
## dbl (2): FIPS_Code, est
## num (2): pop_2010, Area_km2
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#View first few lines of the data frame
head(ny_wiki)
## # A tibble: 6 × 9
## county CountyName FIPS_C…¹ Count…² est forme…³ named…⁴ pop_2…⁵ Area_…⁶
## <chr> <chr> <dbl> <chr> <dbl> <chr> <chr> <dbl> <dbl>
## 1 albany Albany 1 Albany 1683 "One o… "James… 304204 1380
## 2 allegany Allegany 3 Belmont 1806 "Genes… "A var… 48946 2678
## 3 bronx Bronx 5 none 1914 "New Y… "Jonas… 1385108 149
## 4 broome Broome 7 Bingha… 1806 "Tioga… "John … 200600 1852
## 5 cattaraugus Cattaraugus 9 Little… 1808 "Genes… "A wor… 80317 3393
## 6 cayuga Cayuga 11 Auburn 1799 "Onond… "The\x… 80026 2238
## # … with abbreviated variable names ¹FIPS_Code, ²CountySeat, ³formed_from,
## # ⁴named_for, ⁵pop_2010, ⁶Area_km2
Now we will join our ny_wiki information data frame to
our county_ny data frame using a left_join()
from the dplyr package. (Consult your dplyr
cheat sheet for more information on available join functions. You may
find these useful as you continue to work with your project data!) We
will also use the mutate() function to calculate a
people_per_km2 variable to express population density.
# Create joined data frame
ny_info <- left_join(county_ny, ny_wiki, by = "county") %>%
# Calculate population density
mutate(people_per_km = pop_2010 / Area_km2)
# Use str() or head() to inspect your new data frame
Using our new dataframe ny_info, we will plot the map of
New York by referencing our basemap and adding a layer that fills county
boundaries based on our new people_per_km variable.
#Create new map object with layer for pop density
ggplot(ny_info) +
# add sf geometry with fill corresponding to population density
geom_sf(aes(fill = people_per_km), color = "#F5DBCB") +
theme_void()
Hmm…all the counties look pretty much the same. The problem is that the population of New York City is so “yuge” that it makes it hard to see differences among the other counties. We can fix this problem by rescaling the gradient. This is similar to performing a log10 transformation on the y-axis to better observe differences among groups.
#Create new map object with layer for pop density
ggplot(ny_info) +
# add sf geometry with fill corresponding to population density
geom_sf(aes(fill = people_per_km), color = "#F5DBCB") +
theme_void() +
# add layer to density map rescaling the gradient by log10
scale_fill_gradient(trans = "log10")
Challenge 3
Based on what you have learned, create a map for the state of New York that color codes counties based on the year they were established.
Consult Google or your ggplot cheat sheet to learn how to change the color gradient for your map.
Consult Google or your ggplot cheat sheet to learn how to change the title of your legend.
When you are satisfied with your map, useggsave("A10_C3.jpg")to export your map image.
Upload your plot to Challenge 3 on Brightspace.
#Create new map object with layer for time established
ggplot(ny_info) +
# add sf geometry with fill corresponding to year of establishment
geom_sf(aes(fill = est), color = "white") +
# change color of gradient and name of legend
scale_fill_gradient(low = "#F5DBCB", high = "#8B5666", name = "Year Established") +
theme_void()
#Export image to figures folder
ggsave("figures/A10_C3.jpg")
## Saving 7 x 5 in image
Let’s try zooming in on our population density map to just look at
the area around Clinton County. When I view this sf object, it will give
me coordinates for the “Bounding box” which I can use to set limits on
the scale of my graph. I will switch back to
theme_minimal() so that we can see the latitude and
longitude.
# Create data frame with rows of only Clinton county
clinton <- ny_info[ny_info$county == "clinton",]
# Then I will use the information to zoom into my population density map
# Focus on coordinates of the Bounding box
clinton
## Simple feature collection with 1 feature and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.03188 ymin: 44.43861 xmax: -73.34433 ymax: 45.01157
## Geodetic CRS: +proj=longlat +ellps=clrk66 +no_defs +type=crs
## ID state county CountyName FIPS_Code CountySeat est
## 10 new york,clinton new york clinton Clinton 19 Plattsburgh 1788
## formed_from
## 10 Washington County
## named_for
## 10 George Clinton\xa0(1739\x961812), fourth\xa0Vice President of the United States\xa0and first and third\xa0Governor of New York
## pop_2010 Area_km2 geom people_per_km
## 10 82128 2896 MULTIPOLYGON (((-74.03188 4... 28.35912
I’m going to use the boundary of this sf object as a guide to zoom
into my population density map using the function
coord_sf().
ggplot(ny_info) +
geom_sf(aes(fill = people_per_km), color = "#F5DBCB") +
theme_minimal() +
scale_fill_gradient(trans = "log10") +
coord_sf(xlim = c(-74.2, -73.3), ylim = c(44.3, 45.1))
Note that I definitely had to tweak the xlim and ylim a bit to get
the map to look the way I want. Being able to see the axes makes this a
bit easier to do. You can always go back and add
theme_void() to remove the axes once you have adjusted the
limits.
Challenge 4
Based on what you have learned, create a map of population density that zooms into the area around Manhattan and Queens. (Hint: Manhattan is in “New York” County, and Queens is in “Queens” County.)
Usetheme_minimal()to view and adjust your xlim and ylim as needed.
When you are satisfied with your map, useggsave("A10_C4.jpg")to export your map image.
Upload your plot to Challenge 4 on Brightspace.
#Find reasonable max and min starting points
nyc <- ny_info[ny_info$county == "new york"| ny_info$county == "queens",]
nyc
## Simple feature collection with 2 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.02615 ymin: 40.54823 xmax: -73.71102 ymax: 40.84044
## Geodetic CRS: +proj=longlat +ellps=clrk66 +no_defs +type=crs
## ID state county CountyName FIPS_Code CountySeat est
## 31 new york,new york new york new york New York 61 none 1683
## 41 new york,queens new york queens Queens 81 none 1683
## formed_from
## 31 One of 12 original counties created in the\xa0New York colony
## 41 One of 12 original counties created in the\xa0New York colony
## named_for
## 31 King\xa0James II of England\xa0(1633\x961701), who was Duke of York and Albany before he ascended the throne of England, Duke of York being his English title
## 41 Catherine of Braganza\xa0(1638\x961705), Queen of England and wife of King\xa0Charles II of England
## pop_2010 Area_km2 geom people_per_km
## 31 1585873 87 MULTIPOLYGON (((-73.92874 4... 18228.425
## 41 2230722 462 MULTIPOLYGON (((-73.76259 4... 4828.403
#Plot population density
ggplot(ny_info) +
geom_sf(aes(fill = people_per_km), color = "#F5DBCB") +
theme_minimal() +
scale_fill_gradient(trans = "log10") +
coord_sf(xlim = c(-74.05 , -73.5 ), ylim = c(40.5 , 40.9))
ggsave("figures/A10_C4.jpg")
## Saving 7 x 5 in image
You will notice that this map does not look great. The built-in maps
often do not do well around coastlines when you zoom in. For coastline
maps, you would be better served by data resources such as NOAA’s GSHHS:
A Global Self-consistent, Hierarchical, High-resolution Geography
Database.
https://www.ngdc.noaa.gov/mgg/shorelines/gshhs.html
So far, all of our mapping has involved vector graphics, or data
based on lines that connect points in x-y space. Mapping data may also
contain raster data, or data stored as pixels at
various resolutions. These are the same types of data stored in
photographs, and are commonly used in remote sensing and in various
types of aerial imagery. Here we will use the terrainr and
raster packages to plot aerial imagery of the Lake
Champlain Basin and the Altona Flat Rock burn site. This code should
work for any georeferenced aerial imagery stored as a .tif. Note that I
also used the terrainr package to download this aerial
imagery from the USGS database. For more information on how to import
imagery, see the complete introduction to terrainr in
Mahoney (2022).
First, we will use the raster::stack() function to
import the data for the Lake Champlain Basin. Note that this is a very
large geographic area, so the data were imported at a resolution of 100
m. This means that each pixel on the map represents a 100 m x 100 m area
of land. Then, we will convert the pixel data to a data frame object
using the as.data.frame() function, and define the names of
the columns using the setNames() function.
# Import raster data from .tif file
lc_raster <- raster::stack("mapfiles/lakechamp_ortho_USGS.tif")
# Convert raster data
lc_raster_df <- as.data.frame(lc_raster, xy = TRUE)
lc_raster_df <- setNames(lc_raster_df, c("x", "y", "red", "green", "blue"))
Take a closer look at the lc_raster_df data frame. It
contains x and y data, which represent longitudes and latitudes of the
center of each pixel, as well as rgb color values that specify the color
of the pixel. Note that these are stored as numeric values, so you can
analyze these or perform calculations as you would any other numeric
variable. Such analysis is beyond the scope of this class but forms the
basis of most remote sensing. Let’s plot the data as a map layer using
geom_spatial_rgb() function in the terrainr
package. Note that we also add the coord_sf("WGS84")
function from the sf package to specify the projection we
want to use. (This plot contains a lot of information, so it will take a
few seconds to generate!)
ggplot() +
# Add raster layer
geom_spatial_rgb(data = lc_raster_df,
# Required aesthetics r/g/b specify color bands:
aes(x = x, y = y, r = red, g = green, b = blue)) +
# Specify coordinate system
coord_sf(crs = "WGS84")
So that’s a great start! Notice that the area in Canada projects as
missing data, which is plotted as a rather unflattering shade of pink.
Let’s adjust the zoom on the map to exclude Canada by adding a
ylim = argument to our coord_sf() function.
While we’re at it, we’ll add the Lake Champlain shape file back to the
map as an outline and remove the x and y axis labels.
ggplot() +
# Add raster layer
geom_spatial_rgb(data = lc_raster_df,
# Required aesthetics r/g/b specify color bands:
aes(x = x, y = y, r = red, g = green, b = blue)) +
# Add lcsf shape file
geom_sf(data = lcsf, color = "#1282A2", fill = NA) +
# Change longitude limits
coord_sf(crs = 4326, ylim = c(43.3, 44.95)) +
# Remove background shading
theme_minimal() +
# Remove x and y labels
labs(x = "",
y = "")
Scale bars are essential if one of the goals of your
map is to allow a user to determine distances between locations.
Similarly, a north arrow should always be included to
indicate the orientation of the map. The ggspatial package
allows you to add a scale bar and north arrow to any ggplot mapping
object quickly, and it has several tools that allow you to easily
customize the scale and appearance of these objects. Let’s add a scale
bar to our Lake Champlain app using the annotation_scale()
function, and a north arrow using the
annotation_north_arrow() function. Note that both of these
functions will default to placing the element in the bottom left corner
of the map, so we will add the location = "tl" argument to
our north arrow function to specify that we would like it to go in the
top left of the map.
ggplot() +
# Add raster layer
geom_spatial_rgb(data = lc_raster_df,
# Required aesthetics r/g/b specify color bands:
aes(x = x, y = y, r = red, g = green, b = blue)) +
# Add lcsf shape file
geom_sf(data = lcsf, color = "#1282A2", fill = NA) +
# Change longitude limits
coord_sf(crs = 4326, ylim = c(43.3, 44.95)) +
# Remove background shading
theme_minimal() +
# Remove x and y labels
labs(x = "",
y = "") +
# Add a scale bar using annotation_scale()
annotation_scale() +
# Add a north arrow using annotation_north_arrow(location = "tl")
annotation_north_arrow(location = "tl")
Load the help file for the annotation functions to look at the
customization options. (hint: run ?annotation_north_arrow
in the console)
Challenge 5
Based on what you have learned, plot the aerial imagery for the Altona Flat Rock stored in the “flatrock_ortho_USGS.tif” file in the mapfiles folder.
Add the “frsf” shape file to the map to outline the perimeter of the 2018 burn, and change the color of the outline so that it’s easy to see.
Add a scale bar and north arrow to your map.
Note: You may need to adjust limits to y and x axis in thecoord_sf()function for map to look right.
When you are satisfied with your map, useggsave("A10_C5.jpg")to export your map image.
Upload your plot to Challenge 5 on Brightspace.
# Import raster data from .tif file
fr_raster <- raster::stack("mapfiles/flatrock_ortho_USGS.tif")
# Convert raster data
fr_raster_df <- as.data.frame(fr_raster, xy = TRUE)
fr_raster_df <- setNames(fr_raster_df, c("x", "y", "red", "green", "blue"))
ggplot() +
# Add raster layer
geom_spatial_rgb(data = fr_raster_df,
# Required aesthetics r/g/b specify color bands:
aes(x = x, y = y, r = red, g = green, b = blue)) +
# Add frsf shape file
geom_sf(data = frsf, color = "yellow", fill = NA) +
# Change longitude limits
coord_sf(crs = 4326, ylim = c(44.857, 44.88), xlim = c(-73.65, -73.62)) +
# Remove background shading
theme_minimal() +
# Remove x and y labels
labs(x = "",
y = "") +
# Add a scale bar using annotation_scale()
annotation_scale() +
# Add a north arrow using annotation_north_arrow(location = "tr")
annotation_north_arrow(location = "tr")
ggsave("figures/A10_C5.jpg")
## Saving 7 x 5 in image
NEXT WEEK: Group Assignment
Based on what you have learned, create a map for your data project.
Consider adding: points and labels for plot locations, color-coding points for plot locations based on treatments, changing the size/color of points to correspond to variables in your dataset.
ENV 422: Please feel free to exchange ideas and code with your other group members.
When you are satisfied with your map, useggsave()to export your map image and upload it to the Assignment 11 dropbox on Brightspace.
Anderson, Eric C., Kristen C. Ruegg, Tina Cheng, et al. 2017. Chapter 8. Making simple maps with R. Case studies in reproducible research: a spring seminar at UCSC. Online. https://eriqande.github.io/rep-res-eeb-2017/map-making-in-R.html
Dunnington, Dewey. 2021. ggspatial: Spatial data framework for ggplot2. Online. https://cran.r-project.org/web/packages/ggspatial/ggspatial.pdf
Frazier, Melanie. 2020. Overview of Coordinate Reference Systems (CRS) in R. Online. https://www.nceas.ucsb.edu/sites/default/files/2020-04/OverviewCoordinateReferenceSystems.pdf
Lovelace, Robin, Jakub Nowosad, and Jannes Muenchow. 2021. Geocomputation with R. Open Source Textbook. CRC Press. https://bookdown.org/robinlovelace/geocompr/
Mahoney, Michael. 2022. A gentle introduction to terrainr. Online. https://cran.r-project.org/web/packages/terrainr/vignettes/overview.html